subroutine f_function(f,p,pi,rhoi)
implicit none

    real(8),parameter:: gamma=1.4d0
    real(8) f,p,pi,rhoi
    real(8) ci !sound velocity
    real(8) ai
    real(8) bv!between var
    
    ci=dsqrt(gamma*pi/rhoi)
    
    if(p<pi)then
        bv=(gamma-1.d0)/(2.d0*gamma)
        ai=(p/pi)**bv-1.d0
        f=2.d0*ci*ai/(gamma-1.d0)
    endif
    if(p>pi)then
        ai=rhoi*ci
        bv=(gamma+1.d0)/(2.d0*gamma)
        ai=ai*dsqrt(bv*p/pi-bv+1.d0)
        f=(p-pi)/ai
    endif
    if(p==pi)then
        f=0        
    endif
    
    
    
    
endsubroutine